home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / multiroots / gnewton.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-12-09  |  5.2 KB  |  227 lines
  1. /* multiroots/gnewton.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21.  
  22. #include <stddef.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <float.h>
  27.  
  28. #include <gsl/gsl_math.h>
  29. #include <gsl/gsl_errno.h>
  30. #include <gsl/gsl_multiroots.h>
  31. #include <gsl/gsl_linalg.h>
  32.  
  33. #include "enorm.c"
  34.  
  35. /* Simple globally convergent Newton method (rejects uphill steps) */
  36.  
  37. typedef struct
  38.   {
  39.     double phi;
  40.     gsl_vector * x_trial;
  41.     gsl_vector * d;
  42.     gsl_matrix * lu;
  43.     gsl_permutation * permutation;
  44.   }
  45. gnewton_state_t;
  46.  
  47. static int gnewton_alloc (void * vstate, size_t n);
  48. static int gnewton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  49. static int gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
  50. static void gnewton_free (void * vstate);
  51.  
  52. static int
  53. gnewton_alloc (void * vstate, size_t n)
  54. {
  55.   gnewton_state_t * state = (gnewton_state_t *) vstate;
  56.   gsl_vector * d, * x_trial ;
  57.   gsl_permutation * p;
  58.   gsl_matrix * m;
  59.  
  60.   m = gsl_matrix_calloc (n,n);
  61.   
  62.   if (m == 0) 
  63.     {
  64.       GSL_ERROR_VAL ("failed to allocate space for lu", GSL_ENOMEM, 0);
  65.     }
  66.  
  67.   state->lu = m ;
  68.  
  69.   p = gsl_permutation_calloc (n);
  70.  
  71.   if (p == 0)
  72.     {
  73.       gsl_matrix_free(m);
  74.  
  75.       GSL_ERROR_VAL ("failed to allocate space for permutation", GSL_ENOMEM, 0);
  76.     }
  77.  
  78.   state->permutation = p ;
  79.  
  80.   d = gsl_vector_calloc (n);
  81.  
  82.   if (d == 0)
  83.     {
  84.       gsl_permutation_free(p);
  85.       gsl_matrix_free(m);
  86.  
  87.       GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0);
  88.     }
  89.  
  90.   state->d = d;
  91.  
  92.   x_trial = gsl_vector_calloc (n);
  93.  
  94.   if (x_trial == 0)
  95.     {
  96.       gsl_vector_free(d);
  97.       gsl_permutation_free(p);
  98.       gsl_matrix_free(m);
  99.  
  100.       GSL_ERROR_VAL ("failed to allocate space for x_trial", GSL_ENOMEM, 0);
  101.     }
  102.  
  103.   state->x_trial = x_trial;
  104.  
  105.   return GSL_SUCCESS;
  106. }
  107.  
  108.  
  109. static int
  110. gnewton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  111. {
  112.   gnewton_state_t * state = (gnewton_state_t *) vstate;
  113.   size_t i, n = FDF->n ;
  114.  
  115.   GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J);
  116.  
  117.   for (i = 0; i < n; i++)
  118.     {
  119.       gsl_vector_set (dx, i, 0.0);
  120.     }
  121.  
  122.   state->phi = enorm(f);
  123.  
  124.   return GSL_SUCCESS;
  125. }
  126.  
  127. static int
  128. gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx)
  129. {
  130.   gnewton_state_t * state = (gnewton_state_t *) vstate;
  131.   
  132.   int signum ;
  133.   double t, phi0, phi1;
  134.  
  135.   size_t i;
  136.  
  137.   size_t n = fdf->n ;
  138.  
  139.   gsl_matrix_memcpy (state->lu, J);
  140.  
  141.   gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
  142.  
  143.   gsl_linalg_LU_solve (state->lu, state->permutation, f, state->d);
  144.  
  145.   t = 1;
  146.  
  147.   phi0 = state->phi;
  148.  
  149. new_step:
  150.  
  151.   for (i = 0; i < n; i++)
  152.     {
  153.       double di = gsl_vector_get (state->d, i);
  154.       double xi = gsl_vector_get (x, i);
  155.       gsl_vector_set (state->x_trial, i, xi - t*di);
  156.     }
  157.   
  158.   { 
  159.     int status = GSL_MULTIROOT_FN_EVAL_F (fdf, state->x_trial, f);
  160.     
  161.     if (status != GSL_SUCCESS)
  162.       {
  163.         return GSL_EBADFUNC;
  164.       }
  165.   }
  166.   
  167.   phi1 = enorm (f);
  168.  
  169.   if (phi1 > phi0 && t > GSL_DBL_EPSILON)  
  170.     {
  171.       /* full step goes uphill, take a reduced step instead */
  172.  
  173.       double theta = phi1 / phi0;
  174.       double u = (sqrt(1.0 + 6.0 * theta) - 1.0) / (3.0 * theta);
  175.  
  176.       t *= u ;
  177.      
  178.       goto new_step;
  179.     }
  180.  
  181.   /* copy x_trial into x */
  182.  
  183.   gsl_vector_memcpy (x, state->x_trial);
  184.  
  185.   for (i = 0; i < n; i++)
  186.     {
  187.       double di = gsl_vector_get (state->d, i);
  188.       gsl_vector_set (dx, i, -t*di);
  189.     }
  190.  
  191.   { 
  192.     int status = GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J);
  193.     
  194.     if (status != GSL_SUCCESS)
  195.       {
  196.         return GSL_EBADFUNC;
  197.       }
  198.   }
  199.  
  200.   state->phi = phi1;
  201.  
  202.   return GSL_SUCCESS;
  203. }
  204.  
  205.  
  206. static void
  207. gnewton_free (void * vstate)
  208. {
  209.   gnewton_state_t * state = (gnewton_state_t *) vstate;
  210.  
  211.   gsl_vector_free(state->d);
  212.   gsl_vector_free(state->x_trial);
  213.   gsl_matrix_free(state->lu);
  214.   gsl_permutation_free(state->permutation);
  215. }
  216.  
  217.  
  218. static const gsl_multiroot_fdfsolver_type gnewton_type =
  219. {"gnewton",                /* name */
  220.  sizeof (gnewton_state_t),
  221.  &gnewton_alloc,
  222.  &gnewton_set,
  223.  &gnewton_iterate,
  224.  &gnewton_free};
  225.  
  226. const gsl_multiroot_fdfsolver_type  * gsl_multiroot_fdfsolver_gnewton = &gnewton_type;
  227.